Migrate variography/pair-sampling and grouped statistics/cosampling, and convolution from xDEM with modular API - #925

This PR migrates and expands two main groups of features from xDEM:

1. Co-sampling (= sampling of multiple geospatial objects at the same locations) and stratified sampling (=subsampling within groups) paired with grouped statistics (i.e. binning or classification directly on geospatial objects).
2. Pair-sampling (=sampling pairs of points within one object) and variography (= estimation of spatial autocorrelation).
   Note that "grouped" statistics includes "zonal" statistics, for instance when the input is a vector file.

Most of these features are moved or adapted from xDEM code. Notably, all implementations support Dask inputs and laziness. In particular, the variography module implements a specific log-lag pairwise sampling scheme for both regular and irregular data to efficiently estimate a variogram from large rasters or point clouds.

Also, the NaN-robust convolution functions of xDEM are also moved here into filters, to support both Numba/SciPy efficiently and tie to filters where relevant.

Otherwise, most changes are distributed into two modules:

• sampling/ now contains a new pairsampling (for variography) and cosampling/stratified (for grouped stats) with a support (common helper functions) in addition to subsampling.
• stats/ now contains a new variography module, and refactors the previous stats.py module (+ some processing code that non-optimally lived in raster/pointcloud classes themselves) into three submodules: reduction module (that applies the statistical reducers, supporting Dask/MP there directly), selection module (that selects appropriate data using routines from sampling/) and grouping module (that contains all the grouping logic). The stats.py module now only contains parent functions used by the raster/point cloud objects.

To support Multiprocessing for statistics, this module also adds a new multiproc/readers.py module that contains helper to read only specific chunks of rasters/points over multiple passes by storing metadata more practically (while other modules always did a full pass over the entire raster, with map_overlap/block).

I also added Dask support for point cloud reprojection, for both Multiprocessing/Dask, which was fairly straightforward (each chunk can be reprojected independently), but required some logic specific to LAS/LAZ/COPC. I used the opportunity to move functions previous packed in pointcloud/base and pointcloud/pd_accessor into a new pointcloud/dataframe.py which contains low-level helper function often reused (mirroring raster/array.py) and a new pointcloud/transformation.py (containing the logic for chunked reprojection, especially for Multiprocessing).

The API now exposes for both raster and point cloud objects/accessors, added to base.py modules of point cloud and raster:
cosample() and stats(by=) for grouped stats (not passing by= computes global stats on the whole object)
pairsample() and variogram()

Finally, by tying these new functions to xDEM uncertainty/coreg module directly (in a parallel PR), a couple small bugs came up (like floating precision differences for interp_points with Dask/MP), hence the small fixes in other modules with small added tests.

Resolves #895
Resolves #876

Context of previous implementations

Those features have been present in different forms in xDEM for a while (specifically for uncertainty/coregistration), and have long been planned to be migrated here in GeoUtils as generic features that can naturally interface with others.
See discussions: GlacioHack/xdem#588, GlacioHack/xdem#378, GlacioHack/xdem#947

For co-sampling, the code is largely migrated from a cosampling.py module initially drafted in GlacioHack/xdem#759, which aimed to replace the pre-processing functions of coregistration/uncertainty propagation (which always require 2 datasets, point or raster; sampled at the same location).

For variography and pair-sampling, the code is inspired by the old spatialstats of xDEM and migrated from a _metricspace.py module initially drafted in GlacioHack/xdem#759 with aim to improve upon Dask-support and efficient pairwise sampling, built on top of SciKit-GStat. We also add variogram conversion across backends (GPyTorch, GSTools and SciKit-GStat) borrowing logic from code I largely wrote in https://github.com/geo-smart/spacetime-elevation.

For grouped statistics, this PR supersedes #668, #774 and #815 after discussions in #895 and elsewhere. It aims to replace nd_binning (SciPy-based) in xDEM used widely for coregistration and uncertainty propagation. The co-sampling logic is the same as above, and the plotting/management of the grouped statistics is inspired from that of xDEM, but the Dask/MP implementation is new to this PR (and probably its biggest piece).

For filters, this PR moves generic SciPy/Numba NaN-supporting convolution from xDEM as discussed here: as discussed here: GlacioHack/xdem#300.

Implementation details

Variography and pair-sampling

For variography, we create our own "light" Variogram object, that is easier to optimize for efficiency on large datasets and allows to inter-operate with other packages (GSTools, GPyTorch, SciKit-GStat).
This is because there is no single geostatistical package in Python that has it all (the kriging packages have better interpretability with empirical variograms, the GP ones are more computationally efficient for applying kriging but more black-box for kernel estimation/visualization). Additionally, the most modular variography Python package (which is SciKit-GStat) has a practical but heavy Variogram class that holds too much information at once (all pairwise distances). So we need a wrapper for variography to be efficient for Dask/MP on large datasets.

Because variography analyzes PAIRS of observations, it is inherently hard to scale on large data without subsampling: a 10,000 x 10,000 raster has 100 million obs, so roughly 100 M x 100 M / 2 = 5 quadrillion pairs of obs; which is impossible to sample on any hardware (and anyway useless because the spatial correlation is often largely consistent in space).

Thus, in GeoUtils, we focus on adding efficient pair-sampling to naturally make the link to large datasets. The pair sampling is done in log-lag space to adequately sample pairs across all distances (otherwise a pure random sample has low probability of having two neighbouring pixels sampled, for instance), and supports Dask input (but not Multiprocessing, it is a bit too complex for it yet). This expands previous work I did in SciKit-GStat. The full detail of the implementations is a bit complex, and it is available in the classes docstring.

Once this pair sampling is done, our wrapper calls SciKit-GStat for empirical binning of the pairs, and then model fit (the computationally cheap part). Here, we have to use some tricks to fake the presence of some class SciKit-GStat attributes that would normally be too large to fit in memory, and make our pair sampling connect to its underlying MetricSpace class. A bit hacky, but it works! 😄 (And we can adjust more cleanly if a 2.0 comes out in SciKit-GStat at some point, to which I might contribute)

The idea is that the new gu.Variogram can be exported to any desired format by the user (to_gstools(), to_gpytorch(). And, later on, within GeoUtils, it can be passed to any interpolation/reprojection/gridding function to trigger kriging: interp_points(method="kriging", variogram=gu.Variogram(...)).
This is not implemented yet, but it is a simple link to GSTools/GPyTorch, and it will make kriging easy, modular, and relevant to the whole package as a core resampling method which can happen at low level (during CRS reprojection, for instance)!

Grouped stats and co-sampling

For grouped stats, we need to do two things: 1/ compare the inputs on a similar spatial support (the object being analyzed, and the grouping objects) and 2/ apply the statistics and have an implementation that works with chunks.

For 1/, inputs can be rasters or point clouds (continuous binning or categorical grouping), or vectors (zonal = categorical grouping). We have to choose a reference spatial support to compare to, which by default is the input object on which the stats are computed (raster or point cloud), but this can be selected using at=. For a vector input, we first create a categorical masking on the reference.

Then, we re-use logic coded in co-sampling.

When comparing raster and point, we have essentially 4 modes (raster_point_mode):

1. Resample the raster grid at point coordinates, for a continuous variable = Interpolation at point coordinates (interp_points)
2. Resample the raster grid at point coordinates, for an area-averaged variable = Reduce around points (reduce_points)
3. Grid the point cloud at raster coordinates, for a continuous variable = Gridding by interpolation at cell coordinates (e.g. grid(method="idw" or "linear")
4. Grid the point cloud at raster coordinates, for an area-averaged variable = Gridding by reduction of points around cell coordinates (e.g. grid(method="average"))

The co-sampling then returns all datasets sampled at the same location following the above modes, and reference coordinates. Dask support here is directly derived from that in interp_points, rasterize, etc. Nothing new except linking to those existing functions.

From there, a couple routines help perform grouping in coordination with stratified sampling: subsampling within each groups. Those are pretty straightforward by reusing the subsampling module.

For 2/, we need to compute the statistics in chunks, and return them for the whole object.
We can do this exactly with a typical split-aggregate-combine workflow for many estimators (mean, STD, RMSE, etc). But other estimators require the full group loaded at once (median, NMAD, etc).

As this is not rocket science (pretty simple logic), we code our own so that it can also interface with other aspects (subsampling, vector geometry considerations etc) and run both in Dask and Multiprocessing.
In order we do these steps:

1. We define a common spatial support based on inputs and user arguments: either a point support, or raster support.
2. We reproject different-CRS inputs to the CRS (and grid for raster) of the common support using Dask/MP support in reproject(), create_mask(), interp_points() and grid().
3. Optionally, if the user passed only a bin count (e.g., 10) for a grouping variable, we need to quickly inspect per chunk to get the min/max values and define the bin edges (simple min, max chunk accumulation to do out-of-memory).
4. We assign an integer ID to every group (lowest integer-type possible, depending on the number of groups), and write on the common support (out-of-memory, looping over all chunks once),
5. We sample the locations for validity (out-of-memory), and optionally subsample them (globally, or stratified = per group), yielding the final sample of the common support used.
6. We compute statistics (=reduce) all values in the final sample, chunk per chunk. Some functions like "mean", "std", can be derived from merging aggregate scalars (e.g. we sum every pixel value per chunk, then divide later by their count), and thus don't need the full group to be loaded in memory. For statistics like median/NMAD/percentiles, one needs the full group in memory at once (which should still be much smaller than the input raster). This is why we add stratified subsampling: We can still subsample a given group to a max value (e.g., 100,000 samples) to ensure statistics like median/NMAD/etc don't blow up the RAM usage. This step requires its new Dask/MP implementation.

Then, done!

Benchmarking tests show the performance is on-par with Flox, so we have a good implementation 🙂.

Filtering and convolution

I have moved the Numba/SciPy dual-convolution implementation of xDEM in GeoUtils. I used the opportunity to link it to our existing filters, so that it is used consistently. It adds more modularity than the old uniform_filter implementation used for mean_filter before the generic_filter of SciPy came out (e.g., we can pass any shape, including circular, instead of just the square of uniform_filter). I also added a minmax filter in Numba, so that all of our filters consistently support both computational engines: SciPy or Numba! 😄